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We analyze the excision strategy for simulating black holes. The problem is modeled by the 
propagation of quasi-linear waves in a 1-dimensional spatial region with timelike outer boundary, 
spacelike inner boundary and a horizon in between. Proofs of well-posed evolution and boundary 
algorithms for a second differential order treatment of the system are given for the separate pieces 
underlying the finite difference problem. These are implemented in a numerical code which gives 
accurate long term simulations of the quasi-linear excision problem. Excitation of long wavelength 
exponential modes, which are latent in the problem, are suppressed using conservation laws for the 
discretized system. The techniques are designed to apply directly to recent codes for the Einstein 
equations based upon the harmonic formulation. 
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I. INTRODUCTION 



Two computational difficulties in simulating the emission of gravitational waves from black hole interactions are 
the existence of a singularities inside the black holes and a consistent treatment of the outer boundary so that the 
waveform of the gravitational radiation can be accurately extracted. Fortunately, the observable waves emanate only 
from the region outside the black holes and can be simulated, without loss of physical content, in a domain from which 
the singularities have been excised by introducing inner boundaries inside the black holes. However, the treatment 
of this domain by Cauchy evolution has had limited success. We present here a quasi-linear model excision problem 
whose mathematical analysis is simple enough to reveal some underlying computational pitfalls and the means to 
avoid them. In particular, we show that an accurate simulation of this model excision problem can be obtained by 
(i) blending together two evolution algorithms, one of which is well-posed near the outer boundary and the other 
well-posed near the inner boundary, where a superluminal evolution is employed, and (ii) incorporating discrete 
conservation laws which control spurious exponential growth. The results are aimed at numerical formulations of the 
Einstein equation as second order quasi- linear wave equations, such as in recent harmonic approaches Q,!^!^- 

The problem can be illustrated by considering the propagation of scalar waves on the fixed background geometry of 
a spherically symmetric Schwarzschild black hole. A choice of coordinates commonly adopted in numerical work for 
excising the singular region are ingoing Eddington-Finkelstein coordinates = {t, r, 9, 0), in which the Schwarzschild 
spacetime metric gap is given by 

g^pdx^dx'' = -{l-—)dt^ + —dtdr + {l + —)dr^ 



r r r 

,2 /) J i,2 



(dr +sin-^6'dr). (1.1) 



Here, the spacetime is foliated by the spacelike Cauchy hypersurfaces t = const, which are used to carry out the 
numerical evolution. The radial coordinate extends from the location of the singularity at r = to the asymptotic 
region r oo, to which the radiation propagates. The black hole is located at r = 2M. In the excision strategy, the 
simulation is carried out in a region Ri < r < R2, where the inner boundary satisfies < i?i < 2AI and the outer 
boundary (introduced for the purpose of a finite coordinate range) satisfies R2 >> 2M. 
The covariant wave equation for a scalar field u is 



gg'-^dpu) = 0, (1.2) 
where g = det (gap) and g^f^ is the inverse metric. In the above coordinates, it takes the explicit form 

r r J, j,z gjj^z V 

2M 2{r-M) 
= 7r<Jtu ^ OrU — 2 cot Oogu. (1.3) 

Although this wave equation is non-singular and can be reduced to symmetric hyperbolic form in the region 
i?i < r < i?2, it has some notable features which complicate a numerical treatment: 
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The appearance of the mixed dtdr derivative. This is an effect of the Eddington-Finkelstein coordinates in which 
the (r, 9, cj)) — const observers are moving relative to the time shcing. In the standard ADM |J| formahsm, this 
gauge effect is described as a non- vanishing shift, described by the shift vector /3" = (0, (3^, 0, 0) where 

= > 0. (1.4) 

The change in sign of the coefficient of in passing across the black hole. This change in sign does not affect 
the hyperbolicity of the wave equation but it creates a horizon at r = 2M across which waves cannot pass to the 
exterior. For r > 2M, the (r, 9, (p) = const-observers move on timelike worldlines; for r = 2M, these worldlines 
are lightlike; and for r < 2M, they are spacelike. Thus, inside the horizon the dt direction becomes spacelike, 
i.e. an evolution based upon grid points with constant r is superluminal inside the horizon. 

The energy density on the Cauchy hypersurface t = const associated with a timelike vector ^" is £ = £,°'Tapn^ , 
where n^ — — {dat) / \/ — gt^^ {d^t)d^t is the unit timelike normal to the Cauchy hypersurface and 

Tap = {dau)dpu - ^gapg^''{dfj^u)d^u (1.5) 

is the energy momentum tensor of the scalar field u. In this paper we use the conserved energy associated with 
the choice ^" = t", defined by f^da = dt, whose flow leaves the Schwarzschild metric invariant. However, the 
corresponding energy density 



£ = . ^ f (1 + —)idtu)^ + (1 - —)idr.u)^ + \{deu)' + -^^d4,u)A (1.6) 
2 A _|_ 2M V r f ?• sm 9 J 

is positive-definite and provides a norm for the scalar field only for r > 2M, where is timelike. Inside the 
horizon, is spacelike and £ represents momentum rather than energy density. This complicates the use of 
the energy method for establishing a stable algorithm in the interior of the horizon. The non-conserved energy 
associated with the choice ^" = does provide a global norm but we will not pursue that direction here. 
Instead, we use mode analysis to treat the inner region. 

These features are further complicated by the consideration of proper boundary conditions for a numerical treatment. 
They are not peculiar to a second differential order treatment. In particular, all these features are prominent in a 
model excision problem based upon a first order symmetric hyperbolic treatment of the linear wave equation H1.3() 
carried out in Cartesian coordinates They are also prominent in another second order discussion of the excision 
problem , in which the energy norm associated with the choice ^" ~ n" is used to treat the stability of the evolution 
algorithm (|2.8|l which we adopt for the superluminal case. 

We study the numerical simulation of the excision problem by introducing a model 1-dimensional system which 
retains all the geometric features of the full problem. In the domain < x < 1, we consider the quasi-linear wave 
equation 

daiV^g'-^^dp'f) = dt{^g'^dp<i>)+d,{^g^''d0'S>) (1.7) 

= -dt{^dt<i>) + 9*(|a.$) + 9.(|at<i>) + d/-^d,<f) 

= (1.8) 

where the coefficients are smooth and satisfy 5 > and > 5 for x = and < 6 for a; = 1. These conditions 
ensure that the i-foliation is spacelike, that the metric determinant satisfies the Lorentzian requirement g < 0, that 
the outer boundary x = 1 is timelike and that the inner boundary a; = is spacelike. The shift vector has component 
= a. We also require, analogous to l|1.4|l . that a > at a; = 0, which ensures that the inner boundary is oriented 
so that all characteristics leave the domain. There is an event horizon at = 6 across which waves cannot pass in 
the outer direction. For initial data satisfying $ > 0, the Cauchy problem for this equation is well-posed. 

The non-linearity is introduced into 1)1. 8|l in order to model another computational difficulty related to the full 
Einstein equations. When harmonic coordinates are introduced to fix the gauge freedom, the Einstein equations 
reduce to quasi- linear equations for the components of the metric whose principle part is identical to (|1.7(l . In these 
harmonic coordinates, there exists pure gauge solutions (flat metrics) 

ds^ =^-dt^ + dx^)+dy^ + dz^, (1.9) 
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where 



(1.10) 



which describing traveUng waves with_profile / > — 1. The simulation of (|1.9|) is one of the standardized tests proposed 
in the "Apples with Apples" project ItIJS'I, whose purpose is to compare different numerical treatments of the Einstein 
equations. The accuracy of the simulation is complicated by the existence of exponentially growing gauge waves 9] 



which solve the same equations and cannot be damped by dissipative techniques. For A « 0, the Cauchy data for 
p. 11(1 approximate the Cauchy data for the test gauge wave 1)1. 9|l . The excitation of this exponential error mode can 
dominate the simulation after a relatively short time. An accurate long term simulation of p. 9(1 requires controlling 
the excitation of ((1.11(1 by some global method. Excellent long term accuracy has been obtained for the gauge wave 
test (with periodic boundaries) using a code llOl based upon a semi-discrete energy estimate derived at the linearized 
level by the summation by parts method [Til ll2LlT3j | . Comparable results for the gauge wave test can also be obtained 
using analogs of the semi-discrete monopole conservation laws applied in Sec. IIVI to the quasi-linear wave equation 
(fO|) . See Jii for details. 

This quasi-linear wave equation closely models the gauge wave problem. Solutions of 1(1.8(1 can be used to construct 
gauge waves via ((1.9(1 . The traveling waves ((1.10(1 and ((1.12(1 are examples. Thus the accurate long term simulation 
of solutions to ((1.8(1 can be complicated by the existence of neighboring solutions which grow exponentially in time, 
in precise analogy with the gauge wave problem. 

We analyze the quasi-linear excision problem in terms of its elementary pieces. A stable algorithm for an initial- 
boundary value problem must also be stable for the corresponding linearized problem with frozen coefficients. This 
problem is considered in Sec.^ We consider three separate regions: (i) behavior of the interior evolution, (ii) behavior 
at a time-like boundary and (iii) behavior at a spacelike boundary. Treatment of a time-like boundary in general 
relativity introduces yet another computational difficulty related to the constraints. In the harmonic formulation 
of the initial-boundary value problem, each component of the metric satisfies a quasi-linear wave equation whose 
principle part is identical to ((1.2(1 . The only well-posed version of the initial boundary value problem which is known 
at present to preserve the harmonic constraints requires a Dirichlet boundary condition for some metric components 
and a Neumann boundary condition for the others [3] ■ For the purpose of future application of our results to harmonic 
evolution in numerical relativity, we concentrate here on such Dirichlet and Neumann boundary conditions, although 
the treatment extends to the general maximally dissipative boundary condition. 

The spatial discretization of the wave equation in second differential order form has many advantages over first 
order formulations, as discussed in [l4 i Il3 . These advantages have been recognized in recent second order harmonic 
approaches to numerical relativity 0, 0, For the case of zero shift, a comprehensive treatment of the wave 
equation in second order form has been given for both Dirichlet (16j and Neumann boundary conditions |l7l |. A major 
challenge of the excision problem is not just dealing with a shift but with a shift that is superluminal. This introduces 
complications even in the absence of boundaries. Standard implicit [isj and explicit evolution algorithms for 
the second order wave equation with shift are unstable in the superluminal case. We treat this Cauchy problem in 
Sec. Ill Al In Sec. IIIBI we give a finite difference algorithm for the subluminal case which is stable for both Dirichlet 
and Neumann boundary conditions. We give a detailed analysis, since there is little treatment of the second order 
wave equation with shift in the literature. Whereas treatment of the Dirichlet condition carries over unchanged, a 
non-zero shift forces a revision in the treatment of the Neumann condition. The algorithm for the subluminal case 
(the outer algorithm) is well suited for treating the outer boundary of the excision problem but becomes unstable 
inside the horizon, where the shift is superluminal. In Sec. Ill Cl we give an alternative finite difference algorithm (the 
horizon algorithm) which is stable in this region and admits a stable spacelike excision boundary. 

We establish stability by using the method of lines to reduce the wave equation to ordinary differential equations in 
time on a spatial grid and then applying two standard techniques: mode analysis and semi-discrete energy estimates 
supplied by the summation by parts technique |20|. The energy method is important because it extends directly to 
the wave equation with non-constant coefficients in a region with inner and outer boundaries. We use the energy 
associated with the t-direction along which the spatial grid propagates. For the region inside the horizon, where this 
energy does not provide a positive-definite norm, we use mode analysis to establish a stable excision algorithm at the 
inner spacelike boundary. 

In Sec. lIIII we treat the linearized excision problem. This involves blending together the outer and horizon algorithms 
obtained in Sec.^] The additional ingredient is the necessity of non-constant coefficients in order to model a domain 



dsl = ^xi-dt"^ + dx^) + dy'^ + dz^ , 



(1.11) 



where 




(1.12) 
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with spacelike inner boundary and timelike outer boundary. The homogeneous Dirichlet and Neumann boundary 
conditions are constructed so that the energy of the semi-discrete system is conserved when the coefficients of the 
wave operator are time independent. In addition, the algorithm is constructed so that it conserves a monopole moment 
associated with the scalar field in the case of time-dependent coefficients. 

The algorithm is implemented as a numerical code using a fourth order Rungc-Kutta time integration. Two code 
tests of the linear excision problem are presented in Sec. IIII ("I one where the data describe a pulse completely inside 
the horizon; the other where the data consist of an ingoing wave entering the outer boundary. The tests confirm the 
accuracy and stability of the finite difference algorithm. 

In Sec. II VI we treat the excision problem for the quasi-linear wave equation. Although semi-discrete energy conser- 
vation is not considered in the nonlinear case, we construct the algorithm so that the scalar monopole moment remains 
conserved. Studies of the gauge wave problem 1^ show that either energy conservation or monopole conservation are 
sufficient to suppress the exponential mode H1.12|l . This is confirmed in the numerical tests of the quasi- linear excision 
problem presented in Sec. IIV Bl 

In Sec. we discuss the application to the harmonic Einstein equations of the algorithms developed here for the 
model excision problem. 



II. ALGORITHMS FOR THE LINEAR WAVE EQUATION WITH CONSTANT COEFFICIENTS 

We consider the wave equation 

Utt — "^aur^t — {b ~ a^)uj.x = , 5 > (2-1) 

obtained from 1)1. 8|l by freezing the coefRcients and making the linear approximation $ w 1 -I- u. (Where no confusion 
arises, we use the abbreviated notation dtu — Ut, etc.) We treat three separate problems: 

1. The interior Cauchy problem with both signs of (6 — a^). 

2. The initial-boundary value problem in the domain < a; < oo with a timelike boundary at a; = 0, i.e. (5— a^) > 0. 

3. The initial-boundary value problem in the domain < a; < oo and a spacelike boundary at a: = 0, i.e. (6— a^) < 0. 

If we can construct stable difference approximations for these three problems then the combined difference approx- 
imation for the full linear excision problem is also stable and the rate of growth is bounded independently of the step 
size. The reduction of the initial boundary value problem to Cauchy and halfplane problems is discussed in jlflj l and 
in Ch. 12 of ^J. 



A. The Cauchy problem 

The Cauchy problem for (|2.1() is well-posed for every fixed value of a and b > 0. The Fourier mode 
satisfies 

Utt - 2iaLuut - (6 - a^)w^ = 0, (2.2) 

with solution u = e"^* where 

\ = iuj{a±Vb). (2.3) 

Thus 5R(A) = and there are no exponentially growing modes 

On a uniform spatial grid x^, = vh, consider the simplest second-order accurate finite difference approximation to 

W := Utt - 2aDoUt - {b - a^)D+D_u = 0, (2.4) 

where DqUu = {u^+i — u^-i) / {2h), D+Uu = [uu+i — Ui,)/h and D-Uu = {uv — Ui,-i)/h. The discrete Fourier mode 
Uu = ue^"^"^ satisfies 

2ias\n(u)h) ^ Alb - a?) sir? (ujh / 2) ^ „ ,„ ^, 

Utt -Ut + ^ '-^^ '-^ii - 0, 2.5 

h /i^ 



5 



with solution u = e^* where 

A = 



— - ±T\/a^ sin2(w/i) + 4(6 - a^) sm^{u;h/2). (2.6) 



For 6 > it follows that 5R(A) = but for b < there are exponentially growing modes, e.g. 



X^±-Va^ -b. (2.7) 

for the shortest wavelength mode ujh = tt. 

In order to remedy this problem, we consider the alternative second order accurate finite difference approximation 

V := uu - 2aDoUt - {bD+D^ - a^Dl)u = (2.8) 

for which the discrete Fourier mode satisfies 

2iasm(ujh) f ibsin^(ujh/2) ~ sm'^(ujh)\ 
utt ^^^"*+( [2 ^ -]u = 0, (2.9) 

with solution u — c^* where 

A = [asin{ujh) ± 2Vb fi\n{ujh / 2)^ . (2.10) 

Now 5ft(A) = for all values of the coefficients, subject to 6 > 0. 

Unfortunately, because the evolution algorithm (|2.8|) involves a wider stencil than (|2.4|) it is difficult to analyze its 
behavior and design a clean algorithm near a timelike boundary. Accordingly, we will use the outer algorithm (|2.4|l 
near a timelike boundary and the horizon algorithm (|2.8|) in a region extending from the excision boundary to the 
exterior of the horizon. 

B. Timelike boundary problem 

We consider the initial-boundary value problem for (|2.1|l with & — > in the domain < x < 00. Without 
restriction we can assume b — = 1 . Thus we consider the system 

Utt - 2autx - Uxx = 0. (2.11) 
We can derive an energy estimate. For real functions, let 

/•OO 

{u,v) — / uvdx, = (u, u) (2-12) 

Jo 

denote the L2-scalar product and norm. The energy of the system is 

1 1"°° 1 
^=2!, i< + ^l)dx = ^[\\utr + \\uxr)- (2.13) 



Integration by parts gives 



^5t||w(|p = [ut.Uxx) + '2'a{ut,utx) 



We obtain energy conservation 



-\dt\\uxf - Ut{Q,t)ux{Q,t) - aul{Q,t). (2.14) 



Et^\dt{\\utf + \\uxf)^Q (2.15) 



if the boundary condition at a; = is given by the Dirichlet condition 

Mt(0,t) = (2.16) 



or the Neumann condition 

u40,t) + aut{0,t) = 0. (2.17) 

Note that the Neumann condition (|2.17|l involves the normal derivative g'^'^da — dx + adt in the rest frame intrinsic 
to the boundary. 

We now discuss the difference approximation 

■.= u^tt-2aDou^t- D+D^u^^O, ly ^ 0,1,2, ... . (2.18) 

Corresponding to (|2.12() we use the discrete scalar product 

h ^ 

{w,v)h = -wqVo + ^w„v^h. (2.19) 

Using Lemma 11.1.1 of 20], one can easily derive the summation by parts rules 

{w,Dav)h = -{Dow,v)h-^{wi+w^i)vo-^{vi+v-i)wo, (2.20) 

{w,D^v)h = -{D+w,v)h - ^{wivo + wqv^i). (2-21) 



From 

1 



^dt\\ut\\l = iut,utt)h = 2a{ut, DoUt)h + (ut, D^D+u)h. (2.22) 



Using H2.20|l with 'w,v replaced by Ut, we obtain 



2a{ut, DoUt)h = ~^("it + M-it)"ot- (2.23) 

Using H2.21|l with w,v replaced by Ut and D^u, respectively, we obtain 

19 1 

{ut,D_D+u)h = -- — \\D+u\\l~-{uitD+uo + uotD^uo) 

= -~\\D+u\\l~~\D+uof-uotDoUo, (2.24) 

i.e. 

^dt (^\\ut\\l + \\D+u\\l + ^\D+uofh^ = - (^Douo + ^{uu + u_u)) uot- (2.25) 

Thus we obtain discrete energy conservation if we use the boundary condition 

uot = (2.26) 



J\fo:=DoUo + ^iuit + u_u))^0, (2.27) 

corresponding to second order accurate approximations to the Dirichlet condition (|2.1()l) or the Neumann condition 
(ITTTIl . 

The Neumann condition involves the ghost point v — ~1, which we eliminate by means of (|2.18(l to obtain 

-^W(i+No = -^Utm + D+uo + auti. (2.28) 
We use (|2.28|l to update the boundary point via 

utm - ^ (d+uo + auti \ = 0. (2.29) 
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In the continuum problem, a homogeneous Neumann boundary condition also leads to conservation of the monopole 
quantity 

/>oo 

Q — {ut — aux)dx. (2.30) 
This carries over to the conservation of the semi-discrete monopole quantity 

oo ^ 

S = h'^{uty - aDou^) + 2 ("to " aD+uo), (2-31) 
j/=i 

i.e. Qt = 0, when H2.29|) is satisfied. This semi-discrete monopole conservation is extended to the non-linear problem 
in Sec. IIVI where it has proved effective at suppressing the long wavelength slowly growing mode ifTT^ . 



C. Spacelike boundary problem < a; < oo. 

We now treat the initial-boundary value problem for (|2.1|l in the halfspace < a; < oo for the case a? > h with 
a > 0, The condition o? > b implies that the boundary a; = is spacelike and the condition a > implies that the 
boundary is oriented so that the characteristics leave the halfspace. Thus no boundary conditions are necessary in 
the continuum problem. 

Since there is no energy estimate of type (|2.15() , we use mode analysis with frozen coefhcients to formulate a stable 
discretization near the boundary. Consider first the continuum system. We study bounded solutions of the form 

u^e'^ ■ u{x) , 5R(s) > 0. (2.32) 

Introduction of this ansatz into ()2.1|l gives the ordinary differential equation 

s^u — 2asux ~ {b — a'^)uxx — 0. (2.33) 

Thus u has the form 

u = e"*((Tie''i"-f (726^=") (2.34) 
where fj,j are the solutions of the characteristic equation 

- + 2asn+ {b~ a^)fi^ ^0, (2.35) 



as / a^s^ , s^{b-a^) _ as ^ sVb 

^1.2- b_a^^^ {b-a^r^ {b-a^y ~ a^-b a^-b- ^'^^^ 

Since a > Vb, all roots satisfy 3?(/ii,2) > for 5R(s) > and there are no bounded solutions. Thus the Kreiss condition 
(see I23, Ch. 10) is trivially satisfied and the problem is well posed. 
As difference approximation for the halfplane problem we use the horizon evolution algorithm (|2.8|) . i.e. 

(u^+l - Uj._l)t (m^+1 - 2w^ + M^_l) 2('">>+2 - 2u^ +U^_2) r, 

"''^ - 2h ^ + " = ° ('-^^^ 

for ly — 2,S, . . .. We need two extra boundary conditions to determine uq and ui. We use the extrapolation conditions 

h^Dluo = 0, h^Dlm = 0. (2.38) 

Now we study bounded solutions of type 

4 

= e^* ct^k'^j, n{s) > (2.39) 
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where Kj are solutions of the characteristic equation 

s^-asfK- -)-6(k-2 + -) + ^(k- -)2 = 0, s = hs. (2.40) 

Lemma 

(i) (|2.4U|) has no solutions with 

|k| = 1 for3fi(S)>0. (2.41) 

(ii) There are exactly two solutions ki, K2 with 

|kj|<1, j = l,2 for3fi(S)>0. (2.42) 

(iii) For s ^ 0, 



lim «i = -(1 - ^) + , /(I - ^)2 _ 1, lini ^2 = _(l _ ^) - J(i - ^)2 _ 1. (2.43) 

Proof. Part (i) follows from the stability of the Cauchy problem, established in Sec. Ill Al For s real, as s — + cxd the 
solutions of H2.40|l with |k| < 1 converge to the solutions of 

+ — ^ = 0, i.e. Ki,2 = ±-^. 2.44 

Correspondingly, the solutions of (|2.40() with |k| > 1 converge to the solutions of 

(P 2s 
5^ + — ^2 = 0, i.e. K3 4 = ± — . (2.45) 

4 ' ia 

Since the four roots k{s) of H2.40|l are continuous functions of s and |k| ^ 1 for 5R(s) > 0, there are exactly two roots 
Ki_2 with |ki^2| < 1 for all s with 5R(s) > 0, thus establishing (ii). For s — 0, we have 



(k_ 1)2^0 or (k + 1)2 - = _^2(1 - ^)k + 1 = 0. (2.47) 



Therefore, 



2&, /, 26,^ 

'*1.2=- 1-— ±V 1-— 2-1, K3,4 = l. 2.48 

A simple perturbation analysis shows that |k3.4| ~ je^^'^'*! > 1 for \s\ « 1, 5R(s) > 0, i.e. K3_4 correspond to the 
solution of H2.35|l . Since > b, we have |ki.2| = 1 for s ^ and H2.43|l follows. This proves the lemma. 
The lemma shows that as = 0-4 = 0. To account for possible double roots we write (|2.39|l as 

= e"*(criK5' + (J2^^2__!!l) ifKi7^«;2, (2.49) 

which becomes 

Ui, — e**(criK5' + 0'2!^K2~^) if Kl = K2- (2.50) 

Here cri,iT2 are determined by the boundary conditions H2.38|l . i.e. 

ai(^i - 1)=* + {{k2 - If - [k, - If) = , 

K2- Hi 

aiKi{Ki-lf^—^^{K2{K2-lf-Ki{Ki-lf) =0 . (2.51) 
K2 — Ki ^ 
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The determinant of this system is 

Det = (ki - 1)3(k2 - 1)^ (2.52) 

We can now prove: 

Theorem. H2.49|l . H2.51|l has only the trivial solution (ci = (72 = 0) for di{s) > 0. Therefore the difference 
approximation (|2.37() . H2.38I) is stable. 

Proof. By H2.52|l . we have ki = 1 or K2 = 1- Then, by 12.4()|l . s = and, by (|2.43|l . ki ^ 1, k2 1, which is a 
contradiction. Therefore the Kreiss condition is satisfied. In Ch. 12 of [20|, it is proved that the Kreiss condition 
implies stability. 

Remark 1. The Ryabenkii-Godonov condition states that there are no bounded solutions of type H2.49|l for Tl{s) > 0. 
The Kreiss condition requires that there are no bounded solutions for Tl{s) > 0. The Ryabenkii-Godonov condition 
is a necessary condition while the Kreiss condition is a sufhcient condition for stability. (See Ch. 12 of 20].) 

Remark 2. Mode analysis is closely connected with the Laplace transform. Therefore the stability estimates are 
expressed in terms of \\u {-,1)11^^(11. As explained in Ch. 12 of [23|, one can use these estimates to obtain standard 
energy estimates. 

Remark 3. The stability theory in jl^ and in Ch. 12 of is only developed for first order systems. However, it 
generalizes directly to second order equations. Also in it is indicated how to prove that the Principle of Frozen 
Coefficients holds. The approximation with smooth coefficients is stable if the system is strictly hyperbolic and the 
Kreiss condition holds. 



D. Time discretization 



The results in See's (|irB|l and l|!rc|l prove that the outer algorithm gives a stable difference approximation for 
the halfplane problem with timelike boundary and that the horizon algorithm gives a stable difference approximation 
for the halfplane problem with a spacelike boundary. In Sec. (jIII|) . these results will be combined to give a stable 
difference approximation for the strip problem with spacelike inner boundary and timelike outer boundary. 

For the time discretization, we use the method of lines. The spatial discretization reduces the problem to a large 
system of stable ODE's 



Introducing 



we obtain the first order system 



Utt = T^ut + ^Bu. (2.53) 



Ut = iv, (2.54) 



V \ 1 / B A 



u ; h \ I \ u 



V 



(2.55) 



We solve the system numerically using a 4th order Runge-Kutta time integrator. The time step is limited by the CFL 
condition determined from the Cauchy problem. 

III. THE LINEAR EXCISION PROBLEM 

We model a pulse of energy propagating into a horizon using the wave equation 

- d^u + dtiad^u) + d^iadtu) + 9, ( {1 - a^)d^u)] = (3.1) 
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on the interval < x < 1, where we set 

5 1 
a=--x, 0<a;<- 

a = a{x) , -<x<- 

13 , , 

«=2' i^^^l' (3.2) 

with a providing a smooth monotonic transition between the inner and outer regions. The underlying metric is 
Lorentzian and 1)3. l|l is hyperbolic. The inner boundary is spacelike, the outer boundary is timelike and there is an 
event horizon at x = 1/4. 

The solution to H3.1|) in the region 3/4 < x < 1 has the exact form 

2 

u = f{t + -x)+g{t-2x) (3.3) 
o 

consisting of an ingoing wave / and an outgoing wave g. Near the inner boundary all waves are ingoing, as can be 
see by freezing the coefhcient a to obtain the solution 

u = fit+^x)+g(t + 4x), (3.4) 

which is valid in the short wavelength limit. 

Our goal is to simulate a wave / incident on the boundary at a; = 1 which propagates across the horizon, while 
being partially backscattered. We also check that a pulse of compact support inside the horizon does not propagate 
into the exterior region. 

A. Linear algorithm with non-constant coefficients 

First we must modify the evolution-boundary algorithm to account for non-constant coefficients, as required for the 
excision problem. We consider the equation 

— dfU + dt{adxu) -f dx{adtu) + dx{cdxu) — 0, (3.5) 

where a — a{t, x) and c — c{t, x). The energy associated with f^da = dt is 



= i ^ (^^? + cul)dx. 



(3.6) 



Note that for c < 0, i? is not necessarily positive. For the case = Cf = 0, energy conservation holds in the form 

dtE = F\x^i - F\x^o (3.7) 

where the boundary flux is 

F = ut{cux + aut). (3.8) 

With this generalization, the homogeneous Dirichlet condition ()2.16|l remains ut ~ but the homogeneous Neumann 
condition (|2.17l) must be modified to 

cux + aut — 0. (3.9) 

We modify the outer algorithm to 

W := utt - dtiaD^u) - Do{aut) - ^D^ UA+c)D+uj - UA^c)D^u] = 0, (3.10) 



where 



A±U = f!^±l±A_ (3.11) 
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The last two terms in (|3.10|l are equal because of the identity 



D_[iA+f)D+g] ^D+[ iA^f)D^g], 



(3.12) 



but it is advantageous to express in a form which is manifestly reflection invariant. 

Consider now semi-discrete energy conservation for the case at = Ct = 0. We modify the semi-discrete version of 
the energy to 



N-l 



(3.13) 



where 



It follows that 



£ = lu^t+ \{A+c){D+ur + \{A^c){D^ur. 



(3.14) 



£t-UtW = Doiaui)- -\^D+ut)D+{aut)~{D-Ut)D-{aut) 
+ \LtD^[{A+c)D+v^ + '^{D+ut){A+c)D+u 
+ '^utD+({A^c)D^u\ +]^{D^ut){A^c)D^u. 



(3.15) 



Following the procedure leading to H2.25(l . the summation of this expression gives the semi-discrete flux conservation 
law 



dtE 



^MtAT (^Wn + {A+cn)D+un + {A^cn)D^un + (A+aAr)ui(^,+i) + (^_aAr)wf(jv-i) 
iuto ^ - hWo + {A+co)D+Uf) + {A^cq)D^uo + {A+ao)uti + (^_ao)wf(_i)^ . 



We express the Neumann condition H3.9|l in the second order accurate form 

1 



{A^c)D+u + (A_c)D_M + {A+a)T^ut + {A^a)T^ut = 0, 



(3.16) 



(3.17) 



where T±fu — fu±i- At the boundaries, Af'j^ and Aq contain ghost points outside the computational domain, which 
we eliminate via the relations 



-Wn+M'n = ^uuN + iA-CN)D-UN + iA-aN)ut{N-i) 
--^Wo+J^o = --^utto + iA+CQ)D+ua. + {A+ao)uti. 



(3.18) 
(3.19) 



Then, requiring that the wave equation be satisfied at the boundary points, i.e. that Wn = Wq = 0, the conservation 
law ()3.16() implies that the Neumann boundary conditions 



Afjy : = ^uttN + {A-Cn)D^un + iA^aN)ut(N-i) = 
Wq : = -^utta + iA+co)D+uo + {A+ao)uti = 



(3.20) 
(3.21) 



are dissipative. 

Energy conservation holds in the continuum theory only when the coefficients of the wave operator are time 
independent. However, in the continuum problem, conservation of the monopole quantity H2.30|l holds for general 
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time-dependent coefBcients in the case of homogeneous Neumann boundary conditions. This carries over to the 
conservation of the semi-discrete quantity 

N^l ^ 

Q = ft, ^ {utu - aDoUv) + -^{utN - aD^UN + uto - aD+uo) (3.22) 
i/=i 

when the Neumann condition H3.17(l has the correction 

JV — Af' + —atD+D_u = 0, (3.23) 

which preserves second order accuracy. As a result of this correction, H3.20|l and (|3.21|) are modified to 

^Wn + A/'at := ^{uttN - ^inD^un) + {A_cn)D_uj^ + {A_a]s)ut{N-\) = (3.24) 
h h 

--Wo+No := --{utm - atoD+uo) + {A+co)D+uo + (A+ao)wti = 0. (3.25) 
Setting = for 1 < u < N — 1, a straightforward calculation then gives 

Qt = ^Wn + + ^Wo - JVn. (3.26) 
Thus, again requiring that Wn = Wq = 0, the boundary conditions 

A/'at = ^{uttN - atND-UN) + {A^cn)D^un + (A_aAr)ut(Ar_i) = (3.27) 

Wo = ~^{utm - atoD+uo) + {A+co)D+uo + {A+ao)uti = 0, (3.28) 

imply Qt = 0. Thus we have established that Qt = in the general case and £"4 = when the coefficients are frozen 
in time. 

Note that when the coefficients are not frozen in time, an energy estimate still applies provided a,at,c,ct are 
bounded functions, with c > 0. For simplicity, consider the case with periodic boundaries. Then for the continuum 
problem 1)3. 5|l and 1)3. 6|l imply 

Et= I {atutu^ + ^ul)dx. (3.29) 







The inequality 2fg < + gives 



atutu^ < — ^(Wf + cul) < Ki£ (3.30) 



where the constant Ki satisfies Ki > \at\/\/c] and 

^ul<\^{uj + cul)<K,S, (3.31) 

where the constant K2 satisfies K2 > \ct\/c. Integration of (|3.29|) then gives Et < {Ki + K2)E so that 

E < Soe(^i+^^)*, (3.32) 
where Eq is the initial energy. A straightforward calculation based upon the summation 

N 

Et = hJ2{St - utW) (3.33) 
1 

leads to a semi-discrete energy estimate analogous to ()3.32(l . 

Finally, in the inner region, we modify the horizon algorithm (|2.8|l to the non-constant coefficient form 

V := utt - dtiaDou) - Do{aut) - ^D^ (^{A+b)D+u^ - (^{A^b)D^u^ + Do{a^DQu) = 0, (3.34) 

where c = b — a? . 
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B. Blending the outer and horizon algorithms 



For the purpose of excision, we need a prescription for switching from the outer algorithm W ^ to the horizon 
algorithm V — 0. We introduce a smooth, monotonic blending function f{c)=f{b—a^), satisfying / = 1 for c < 
and / = for c > 1/2. Referring to (|3.1U|) and (|3.34|l . we then use the blended evolution algorithm B = 0, where 

B = W - D^(^{A+{fa^))D+u^ + Do(^fa^DQU^ ^ W + ^D+D^{fa^D+D^u), (3.35) 

which reduces to W near the outer boundary and to V inside the horizon. 

When dtg"^^ = 0, the blended algorithm satisfies a semi-discrete energy conservation law. For periodic boundary 
conditions (identifying the points ly ^ and i/ = N), this takes the form dtEs = where 

N 

EB = hY,^B, (3.36) 
1 

with 

8B=£+^-^{D+D^uf. (3.37) 

o 

(Here, as before, £ is not necessarily positive inside the horizon.) 

For the boundary conditions of the excision problem, the energy (|3.13|) is modified in the same way, 

EB^E + h^^-^{D+D^uf. (3.38) 

Since the support of the blending function / is isolated from the outer timelike boundary, the semi-discrete energy 
fiux through the outer boundary remains given by the first term in H3.16|l . 



C. Tests of the linear excision algorithm 

In order to validate the linear algorithm a set of test-runs was performed where the outer boundary condition at 
X — 1 was either Dirichlet or Neumann, while at x = the extrapolation conditions (|2.38(l were applied. The gridstep 
was chosen iohe h — l/{p ■ 2000), p = 1,2,4 with the time-step set to At — h/lQ. The background coefficients were 
set according to with a{x) = 245/4 - 513 a; -I- 1728x2 - 2880^3 -l- 2368 x** - 768 x^ so that differentiability 
results. 

In the first set of runs the initial data &i t = were a pulse of compact support inside the horizon with 
Mt(0, x) = and 

u(0,x) = 0.5 X [4^(1- O]'^ for xq < x < xq 0.2 (3.39) 
u(0, x) = elsewhere, 

with ^(A) = (A — 0.025)/0.2 . Homogeneous Dirichlet data w(i, 1) = was given at the outer boundary. At the initial 
location of the pulse inside the horizon, both characteristics point to the left. The condition ut(0, x) = implies that 
the initial data is the superposition of two pulses, whose amplitudes have opposite signs, which propagate along these 
two characteristics. This is clearly seen in Fig. ^ where at first the "faster" mode (with small negative amplitude) 
falls into the excised region, while the "slower" mode (with positive amplitude) propagates in the same direction. We 
also monitor the convergence factor for the error modes that have propagated outside the horizon, defined in terms 
of the L2 norm over the interval xr — 0.25 < x < 1 as 

C(x > X,) = log, f . (3.40) 

Second order accuracy would imply C{x > xh) > 2 since the analytic solution vanishes outside the horizon. 

In the region outside the horizon, the simulation consists purely of short wavelength error. At first there is a high 
convergence rate but after about 6 crossing times the short wavelength error is unresolved on the grid. The addition 
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FIG. 1: Sequence of snapshots showing the initial pulse falling into the excised region. The left graph shows the early behavior 
as the faster mode (with negative amplitude) propagates through the excision boundary. The right graph shows the slower 
mode making its way through the excision boundary, with no observable trace remaining at t = 2.0. 
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FIG. 2: Plot of the convergence factor C{x > xh) as a function of time for the test run with an initial pulse inside the horizon. 
It is evident that the addition of dissipation is necessary for long term convergence. 



of dissipation to the evolution scheme greatly extends the time for which convergence holds in the outer region, as 
shown in Fig. El Dissipation is added by modifying (|2.55() according to 



ut-^ut + (uD+D^FD+D^u, 



(3.41) 
(3.42) 



where we set e„ = e„ = 0.1 • /i^ and F is a smooth function with F — hva the interior and F = near both boundaries. 

In the second set of runs, the initial Cauchy data were set to zero, u(0, x) = Mt(0, x) = 0, and a wave was introduced 
through the outer boundary x — Ihy prescribing inhomogeneous Neumann data 



3 1 



(3.43) 



Figure shows snapshots of the evolution for a wave generated by boundary data q{t) consisting of a single pulse. 
The pulse enters the outer boundary in the incoming mode of 1)3. 3|l with characteristic velocity —3/2, propagates 
across the blending region and leaves the grid at the excision boundary. The figure also shows the remnant signal 
when no dissipation is added to the code. After the initial pulse has entered the outer boundary, the Neumann data is 
homogeneous and reflects any signal propagating to the right but signals propagating to the left can leave the system 
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FIG. 3: Sequence of snapshots showing a pulse propagating across the grid. The left graph shows the propagation of the main 
signal, while the right graph illustrates the behavior of the residual error. The short wavelength error can be eliminated by 
dissipation. 



across the excision boundary. At t — 1, short wavelength noise fills the region inside the horizon. By t — 5, the short 
wavelength error has either propagated through the excision boundary or has been converted into longer wavelengths. 
At t = 10, no visible signal remains. The short wavelength noise aX t = 1 can be effectively eliminated by a small 
amount of dissipation. 

We checked long term global convergence by prescribing a periodic wave entering though the outer boundary with 
the inhomogeneous Neumann data 



q{t) = 0, for t < 1/150 



q{t) = (cd^ + adt^ Asin'*(7r[x-1.01 + 1.5t]) 



,for i> 1/150. 



(3.44) 



Again the initial Cauchy data was set to zero. The C"^ boundary data 13.4411 is initially set to zero to provide C° 
consistency with the Cauchy data. No artificial dissipation was added in this test. 

We measured Cauchy convergence of the numerical solution by monitoring the convergence factor 



C = l0g2 



(3.45) 



\Up=2 - Up=4||2, 

For the given data and grid sizes, the code displays second order accuracy to within less than one percent. Figure 0] 
plots the convergence factor as well as a snapshot of the residual finite difference error at the end of the simulation, 
at t = 100. 



IV. THE QUASI-LINEAR EXCISION PROBLEM 

A. Algorithms for the quasi-hnear wave equation 

We extend our algorithms to the quasi-linear wave equation H1.8|l in such a way that the outer algorithm continues 
to obey a semi-discrete monopole conservation law in the case of homogeneous Neumann boundary conditions. We 
again set c = b — a? and treat the case of time-dependent, non-constant coefficients. We generalize the outer algorithm 
(|3.10() to the finite-difference form 

^^=^*(^**) -^^-((^+|)^+'^) -^^+((^-1)^-*^) ^t^' (4.1) 

and generalize the horizon algorithm H3.34|l to 

V d, (^U^ -a, (^^i^o*) -D, (^1$,) - (^(A+^)d^$) - \d+ (^(A_ A)z?_cE>^ +D, (^^?o$) - 0. (4.2) 
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FIG. 4: Convergence plots for the test run with a periodic wave introduced through the outer (timehke) boundary. On the 
left, the plot of the convergence factor C vs. time indicates good second order convergence up to the end of the simulation 
at t — 100. The quality of the long term performance is reinforced by the right graph which shows that, at the end of the 
simulation, the rescaled error profiles are in nearly perfect agreement. No artificial dissipation was introduced. 



We blend these algorithms using the quasi- linear version of H3.35|l . i.e i? = with 

B^W-D^ (^{A+^)D+^^ + Do (^Do^^ = W+ ^D+D^ (^l^D+D^^ ) . (4.3) 

As the boundary condition for the horizon algorithm at the spacelike excision boundary, we use the extrapolation 
conditions H2.38|l . now written as 

h^D\<^o = 0, h^D\^i = 0. (4.4) 

Similarly, as homogeneous Dirichlet condition for the outer algorithm at a timelike outer boundary we retain the linear 
form $t = 0. We formulate the homogeneous Neumann condition for the outer algorithm to establish semi-discrete 
monopole conservation corresponding to the continuum conservation law 

-^'"-"^ («) 

We proceed as follows. 

For periodic boundary conditions, (|4.1|l implies the semi-discrete conservation law dtQ — Q where 

^-%{^-'-^} («) 

This conservation law can be extended to a homogeneous Neumann boundary condition by adding boundary terms 
to H4.6|l in the form 

2 V ^0 *o )' 

We generalize the Neumann condition H3.23|l to the non-linear form 
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FIG. 5: Sequence of snapshots showing the initial pulse falling into the excised region, evolved with the quasi-linear system. 
The left graph shows the faster mode (with negative amplitude), while the right graph shows the slower mode making its way 
through the excision boundary. Grid-resolution for this simulation was h = 1/8000. 

As a result of the nonlinearities, (|3.24|) and (|3.25|) are modified to 

-Vo+AAo: = -'l(dA-dAD+^o)+{A+^)D+^o + {A+^)<Pn. (4.10) 

Setting Wi, = for l<i'<N — 1, sl straightforward calculation leads to the nonlinear version of H3.26|l . i.e. 

Qt = ^Wn + Mn + ^Wo - Mo. (4.11) 
We implement the homogeneous Neumann condition in the form 

Then the requirement that Wn = Wo = leads to the semi-discrete monopole conservation Qt — 0- 



B. Tests of the quasi-linear excision algorithm 

The excision tests performed vi^ith the quasi-linear code were based upon the ones performed with the linear code, 
with the (time independent) coefficients of the wave operator having spatial dependence determined by H3.2I) and with 
the same Cauchy data or boundary data now prescribed for ($ — 1) as were prescribed for the linearized field u in 

sec. irrni 

The initial data for the first set of runs consisted of a pulse of compact support inside the horizon with homogeneous 
Dirichlet data at the outer boundary. The initial peak amplitude of the pulse (|3.39|) was ^max = 1-5, putting it in 
the nonlinear regime. Since the characteristic speeds determined by the background metric are the same as in the 
linearized case, the snapshots of the quasi-linear evolution shown in Fig. are qualitatively similar to the linearized 
ones in Fig. ^ There is no visible propagation of the signal into the region outside the horizon. The convergence 
factor C{x > xh) for the error modes that have propagated outside the horizon is defined by replacing m by ($ — 1) 
in (|3.40l) . As in the linear case, the short wavelength error in the exterior region is unresolved on the grid. The plots 
in Fig. show that dissipation must be added to the code to obtain long term convergence. 

In the second set of non-linear runs, the initial Cauchy data were set to zero and a wave introduced through the 
outer boundary by prescribing inhomogeneous Neumann data. Figure [T] shows snapshots of the evolution for a wave 
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FIG. 6: Plot of the convergence factor C{x > xh) vs. time, for the non-linear test run for an initial pulse inside the horizon. 
Addition of dissipation is necessary for long term convergence. 
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FIG. 7: Sequence of snapshots showing a pulse propagating across the grid with resolution h — 1/2000. The left graph shows 
the propagation of the signal itself, while the right graph illustrates the behavior of the residual modes. The long wavelength 
mode, which is apparent ai t = 1, decays by f = 5 and no visible signal remains at t = 10. 



consisting of a single pulse entering the outer boundary. No artificial dissipation was used. The main signal behaves 
qualitatively similar to the linearized case in Figure 13 but the error modes behave differently. At t = 1 the short 
wavelength error again fills the region inside the horizon but now long wavelength error has also been excited which 
extends to the outer boundary. This corresponds to mixing in the exponential modes 

q, = eHt-^)^ (4.14) 

which are exact solutions (for any A) to the quasi-linear wave equation H1.8(l in the region x > 3/4. These modes are 
consistent with the homogeneous Neumann boundary condition (|3.43|) in effect after the pulse has entered the system 
and q{t) vanishes. However, the mode (|4.14() is not excited with positive A and decays by f = 5. No visible signal 
remains at t = 10. Thus the discrete monopole conservation built into the quasi-linear system suppress exponentially 
growing long wavelength modes. 

A further test of long term performance is the introduction of a periodic wave through the outer boundary by pre- 
scribing the inhomogeneous Neumann data (|3.44l) . We measured Cauchy convergence by monitoring the convergence 
factor C (|3.45|) applied to the numerical solution for <i>. Figure |H1 plots C as well as a snapshot of the residual finite 
difference error at t = 100, when the simulation was ended. The code displays second order accuracy to better than 
1 percent. No artificial dissipation was used in the test. 
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FIG. 8: Convergence plots for the non-linear test run with a periodic pulse introduced through the outer (timelike) boundary. 
On the left, the time dependence of the convergence factor C indicates good second order convergence. This result is reinforced 
by the right graph where, at the end of the simulation nt t = 100, the rescaled error profiles show a nearly perfect agreement. 
No artificial dissipation was used in this test. 



We have based successful simulations of the model excision problem for a 1-dimensional quasi-linear wave equation 
on evolution and boundary algorithms which were proved to be stable in the linear case. The quasi-linear algorithm 
incorporates a conserved scalar monopole moment which suppresses the excitation of long wavelength exponentially 
growing modes, such as (|4.14|l . which are latent in the system when Neumann boundary conditions are used. The 
global excision strategy involves matching an exterior evolution algorithm to a horizon algorithm for evolution of 
the interior region. The exterior algorithm admits simple implementation of an outer boundary condition and it 
obeys a discrete version of flux-energy conservation which guarantees stability. Inside the horizon, this energy is not 
positive-definite and the exterior algorithm is unstable. The horizon algorithm is stable in the interior region and 
admits a stable extrapolation algorithm to update the spacelike excision boundary. The choice of horizon algorithm 
presented here is based upon centered differencing. Other choices are possible in the model 1-dimensional problem. In 
particular, a horizon algorithm based completely upon one-sided differencing is stable and requires no extrapolation 
at the inner boundary. This would be a very attractive feature if it could be extended to the higher dimensional 
problem. 

Preliminary investigations indicate that the computational algorithms presented here can be taken over in a fairly 
straightforward manner to the harmonic formulation of the full 3-dimensional gravitational problem. In that case the 
Einstein equations reduce to to the quasi-linear form 



where S''"' consists of lower differential order terms which do not contribute to the principle part. These equations 
have a well-posed initial-boundary value problem when the boundary data is explicitly prescribed 2]. 

The principle part of (|5.1(l can be treated in the same way as a quasi-linear wave equation (|1.2() for each component 
of the metric. From the point of view of designing a code, the two main differences from our model excision problem are 
(i) the the 3-dimensional nature of the full problem and (ii) the constraints which enter into the boundary conditions. 

The evolution algorithms readily extend to 3-dimensions. When S*^" is neglected, the flux conservative form of 
(|5.1() can be carried over to the 3 dimensional discretized system to obtain discrete analogues of the scalar wave 
monopole conservation. There results a conserved quantity representing a spatially averaged rate of growth for each 
metric component, which has monopole symmetry for the g** component, dipole symmetry for the g*' components 
and quadrupole symmetry for the g^^ components. The exact form of the conserved quantities depends upon how 
the principle part is split off in conservative form, with H5.1|) just one possibility. See the discussion in for how 
a splitting analogous to the logarithmic form (|1.8|l for the quasi-linear scalar problem leads to accurate long term 
simulation of the gauge wave problem; it remains to be seen if this splitting is effective in more general spacetimes. 
When S^'^ is neglected and the metric coefficients in the wave operator are frozen, semi-discrete energy conservation 
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(5.1) 
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also results for each metric component in the case of periodic boundary conditions. With a cubic grid boundary, the 
results for the quasi-linear wave equation generalize to dissipative Dirichlet and Neumann boundary conditions for the 
metric components. The chief complication is the application of Neumann conditions at the edges and corners of the 
cube. Here among possible approaches are the use of the summation by parts technique, as in 12], or the introduction 
of a smooth boundary which is treated by interpolation, as in 16, .17j. 

The constraints on the boundary data present another complication. The constraints on the system (|5.1|l . which 
guarantee that the solutions satisfy Einstein's equations, are the harmonic coordinate conditions 



where are explicit harmonic driving functions. These constraints are satisfied by the solutions of (|5.1|) if a certain 
mixture of Dirichlet and Neumann boundary conditions are used for the metric components 2J. In that case, the 
constraints are satisfied when homogeneous boundary data are given and the initial-boundary problem for Einstein's 
equations is well-posed. In the case of inhomogeneous constraint-preserving boundary data it is not known whether 
the system remains well-posed. 

It should be emphasized that, in addition to the computational difficulty, there are analytic and geometric problems 
in treating black holes. There is the possibility of exponentially growing perturbations to the analytic problem in 
the inner region near the excision boundary. There is the nonlocal nature of the true (null) event horizon. Although 
the spacelike apparent horizon can be traced out as the evolution proceeds, it is impossible to locate the null event 
horizon until the exterior evolution is complete. The matching of an inner horizon algorithm to an outer algorithm 
must be carried out across an artificial horizon which results from a given choice of shift, as in our model problem. In 
the black hole excision problem, such an artificial horizon could be defined in terms of the hypersurface across which 
the evolutiorigoes from superluminal to subluminal. 

Pretorius jJl has recently obtained promising results using a second order harmonic code to simulate a black hole 
by means of excision. We are now in the process of applying the new techniques presented here to this problem. 
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